La classificazione dei movimenti sismici è un problema di grande rilevanza per la sismologia e la protezione civile. Non tutti i segnali sismici registrati sono causati da terremoti: eventi come esplosioni artificiali, frane sottomarine o tsunami generano onde sismiche che possono essere erroneamente interpretate. Distinguere in modo automatico e affidabile tra un terremoto e altri tipi di eventi è fondamentale per migliorare la risposta ai disastri, ridurre i falsi allarmi e affinare i modelli di monitoraggio sismico.
Questo progetto propone vari modelli di classificazione in grado di riconoscere i terremoti da altre sorgenti sismiche, sfruttando tecniche di machine learning applicate a dati sismici. L’obiettivo è migliorare l’accuratezza dell’analisi sismologica, supportare le decisioni delle autorità competenti e contribuire a una gestione più efficace delle emergenze.
library(tidyverse)
library(heatmaply)
library(mice)
library(visdat)
library(caret)
library(gridExtra)
library(MASS)
library(tree)
library(class)
library(plotly)
library(pROC)
library(sf)
library(terra)
library(leaflet)
library(rnaturalearth)
library(rnaturalearthdata)
library(UpSetR)
library(naniar)Importiamo vari dataset scaricati da:
https://earthquake.usgs.gov/earthquakes/search/ , il sito è una piattaforma ufficiale del United States Geological Survey (USGS): l’ente scientifico del governo degli Stati Uniti responsabile del monitoraggio e della ricerca sui diversi movimenti sismici. Dal portale è possibile estrarre dati su movimenti del suolo in base al periodo selezionato. Scarichiamo diversi dataset per avere maggiori informazioni e allenare meglio i nostri classificatori.
df1 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (1).csv")
df2 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (2).csv")
df3 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (3).csv")
df4 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (4).csv")
df5 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (5).csv")
df6 <- read.csv("C:\\Users\\Tommaso\\OneDrive\\Desktop\\query (6).csv")
df <- rbind(df1,df2,df3, df4, df5, df6)
df <- unique(df)
head(df)Il dataset si presenta con 22 variabili e 74602 osservazioni.
latitude: La latitudine dell’epicentro, in gradi (valori negativi indicano latitudini sud).
longitude: La longitudine dell’epicentro, in gradi (valori negativi indicano longitudini ovest).
depth: La profondità dell’evento in chilometri.
mag: La magnitudo stimata dell’evento, che misura l’energia rilasciata.
magType: Il tipo di magnitudo calcolata (ad es. Mw, Ml, Ms, ecc.) che indica il metodo o l’algoritmo utilizzato.
nst: Il numero di stazioni sismiche utilizzate per calcolare la magnitudo o la localizzazione dell’evento.
gap: Il “gap azimutale” in gradi, ovvero il più grande angolo tra stazioni sismiche contigue, che fornisce un’indicazione sulla copertura sismica attorno all’epicentro.
dmin: La distanza orizzontale (in gradi) tra l’epicentro e la stazione sismica più vicina; un valore minore tende a indicare una migliore localizzazione.
rms: Il valore “root mean square” (RMS) dei residui dei tempi di arrivo, espresso in secondi; serve a valutare la bontà della soluzione di localizzazione.
net: Il codice della rete sismica (ad esempio, “us”, “ci”, “ak”, ecc.) che ha fornito i dati per l’evento.
id: Un identificatore univoco dell’evento assegnato dalla rete o dal sistema USGS.
updated: Il timestamp (in millisecondi) dell’ultima modifica o aggiornamento dei dati relativi all’evento.
place: Una descrizione testuale della località, spesso indicante la città o la regione più vicina.
type: Il tipo di evento sismico (ad es. “earthquake”, “quarry”, ecc.).
locationSource: L’origine del dato di localizzazione (la rete o il centro che ha determinato la posizione).
magSource: L’origine dei dati relativi alla magnitudo.
horizontalError: L’errore stimato della posizione orizzontale, espresso in chilometri.
depthError: L’errore stimato sulla profondità, in chilometri.
magError: L’errore (stima della deviazione standard) associato alla magnitudo.
magNst: Il numero di stazioni specificamente utilizzate nel calcolo della magnitudo.
status: Lo stato dell’evento, ad esempio “automatic” se segnalato automaticamente o “reviewed” se successivamente verificato manualmente.
Possiamo notare nella correlation heatmap che le correlazioni non sono in media molto alte: infatti non superano mai lo 0,6.
#Selezioniamo Variabili Numeriche
dfnum <- df |> dplyr::select(where(is.numeric))
#Una prima cor matrix senza NA
R <- cor(na.omit(dfnum))
heatmaply_cor(
R, cellnote = R, dendrogram = "none", cellnote_size = 10,
cellnote_textposition = "middle center",
main = "Correlation Heatmap:",
colors = colorRampPalette(c("#8C6D31", "white", "#31A354"))(200)
)Filtriamo solo le variabili numeriche significative per la nostra analisi.
Osserviamo come si distribuisce la variabile type.
ggplot(data = earth |> dplyr::filter(type != "earthquake")) +
geom_bar(aes(x = type), fill = "#A1D99B") + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + coord_flip()Generiamo la variabile target oggetto di analisi: il nostro obiettivo è di svolgere un’analisi di classificazione riguardante un movimento sismico e capire se esso è stato causato da un terremoto oppure da un altro evento, come un’esplosione. Dividiamo quindi la variabile target in 2 classi: “earthquake” e “non-earthquake”.
earth$type <- as.factor(ifelse(earth$type == "earthquake","earthquake", "non-earthquake"))
earth_noNA$type <- as.factor(ifelse(earth_noNA$type == "earthquake","earthquake", "non-earthquake"))
print(earth |> group_by(type) |> summarize(n()))## # A tibble: 2 × 2
## type `n()`
## <fct> <int>
## 1 earthquake 61497
## 2 non-earthquake 13105
Esaminiamo la distribuzione della variabile type appena creata.
Utilizzando la latitudine e la longitudine rappresentiamo su una mappa di mercatore interattiva dove sono situati alcuni movimenti sismici presenti nel nostro dataset. Esaminando la localizzazione delle singole osservazioni, possiamo notare come gran parte dei dati classificati come “non-earthquake” hanno una longitudine e una latitudine molto simile (si trovano tutti a nord-ovest degli Stati Uniti), dunque escludiamo a priori le due variabili dai nostri modelli per evitare che la classificazione venga eccessivamente influenzata dalla posizione geografica (visualizziamo solo un campione dei nostri dati per facilitarne la comprensione).
#preparazione al grafico
sample <- earth_noNA |> filter(type == "non-earthquake")
sample1 <- earth_noNA |> filter(type == "earthquake")
plot <- rbind(sample[1:500,], sample1[1:500,])
get_marker_color <- function(type) {
color <- ifelse(type == "earthquake", "red", "green")
return(color)
}
leaflet() |>
addTiles() |>
setView(lng = 0, lat = 0, zoom = 2) |>
addMarkers(lng = plot$longitude , lat = plot$latitude,
popup = paste("Magnitudo: ", plot$mag),
icon = icons(
iconUrl =
paste0("https://raw.githubusercontent.com/pointhi/leaflet-color-markers/master/img/marker-icon-", get_marker_color(as.factor(plot$type)), ".png"), iconWidth = 25, iconHeight = 41,
iconAnchorX = 12, iconAnchorY = 41,
popupAnchorX = 0, popupAnchorY = -41
)) Rappresentiamo nel seguente scatterplot l’andamento di nst in funzione di mag, ovvero l’andamento del numero delle stazioni di rilevamento di movimenti sismici in funzione della magnitudo. Si può vedere che con l’aumentare di mag aumentano anche il numero di stazioni che rilevano i movimenti.
ggplot(data = earth_noNA, aes(x = mag, y = nst)) +
geom_point(col = "#A1D99B", size = 2, alpha = 0.5) + geom_smooth(col = "darkgreen")## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
b1 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = depthError), fill = "#31A354", alpha = 0.6) +
labs(title = "depthError", y = "")
b2 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = horizontalError), fill = "#31A354", alpha = 0.6) +
labs(title = "horizontalError", y = "")
b3 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = depth), fill = "#31A354", alpha = 0.6) +
labs(title = "Depth", y = "")
b4 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = nst), fill = "#31A354", alpha = 0.6) +
labs(title = "nst", y = "")
b5 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = gap), fill = "#31A354", alpha = 0.6) +
labs(title = "gap", y = "")
b6 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = dmin), fill = "#31A354", alpha = 0.6) +
labs(title = "dmin", y = "")
b7 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = rms), fill = "#31A354", alpha = 0.6) +
labs(title = "rms", y = "")
b8 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = mag), fill = "#31A354", alpha = 0.6) +
labs(title = "Magnitude", y = "")
b9 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(y = magError), fill = "#31A354", alpha = 0.6) +
labs(title = "magError", y = "")
grid.arrange(b1,b2,b3,b4, ncol = 2)Il grafico mostra i boxplot delle nostre variabili di interesse: è utile ai fini delle nostre analisi successive notare che alcune di esse hanno outlier molto evidenti.
#Istogrammi
b1 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = latitude, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "Latitude", y = "")
b2 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = longitude, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "Longitude", y = "")
b3 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = depth, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "Depth", y = "")
b4 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = nst, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "nst", y = "")
b5 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = gap, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "gap", y = "")
b6 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = dmin, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "dmin", y = "")
b7 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = rms, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "rms", y = "")
b8 <- ggplot(data = data.frame(earth_noNA)) +
geom_histogram(aes(x = mag, y = after_stat(density)), fill = "#A1D99B", alpha = 0.6, col = "#8C6D31") +
labs(title = "Magnitude", y = "")
grid.arrange(b1,b2,b3,b4,b5,b6,b7,b8, nrow = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Risulta particolarmente evidente che nessuna delle distribuzioni sembra seguire una andamento normale (vengono subito a mancare le ipotesi per applicare un modello discriminante lineare e anche per QDA).
#Boxplot per classi
b3 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = depth, fill = type), alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "Depth", y = "") + theme(legend.position = "none") +
coord_cartesian(ylim = c(0, 110))
b4 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = nst, fill = type), alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "nst", y = "") + theme(legend.position = "none")+
coord_cartesian(ylim = c(0, 110))
b5 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = gap, fill = type),, alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "gap", y = "") + theme(legend.position = "none")
b6 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = dmin, fill = type), alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "dmin", y = "") + theme(legend.position = "none")+
coord_cartesian(ylim = c(0, 15))
b7 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = rms, fill = type), alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "rms", y = "") + theme(legend.position = "none")+
coord_cartesian(ylim = c(0, 3))
b8 <- ggplot(data = data.frame(earth_noNA)) +
geom_boxplot(aes(x = type,y = mag, fill = type), alpha = 0.6) +
scale_fill_manual(values = c("#8C6D31", "#31A354"))+
labs(title = "magnitude", y = "") + theme(legend.position = "none")
grid.arrange(b3,b4,b5,b6, ncol = 2)Le classi sembrano mostrare varianze differenti in molte delle variabili di interesse.
Esaminiamo gli outlier più evidenti dai boxplot non condizionati:
## latitude longitude depth mag nst gap dmin rms magNst horizontalError
## 1 43.66650 -121.3963 3.301 2.6 9 139 0.1901 0.92 1 45.637
## 2 46.70500 -122.7927 7.576 3.1 8 114 0.3582 0.39 4 45.660
## 3 43.90833 -120.9353 6.153 3.1 8 133 0.1073 0.45 3 61.651
## 4 46.70150 -122.7680 0.135 3.0 11 84 0.2121 0.15 4 1.494
## 5 46.69267 -122.7748 0.212 3.0 11 150 0.2123 0.14 3 0.756
## 6 46.32167 -117.5747 7.849 2.7 15 243 0.6691 0.42 5 76.120
## 7 44.37833 -123.3005 3.699 2.5 17 139 0.5222 0.23 7 1.245
## 8 47.70600 -119.8953 17.753 3.1 6 339 0.8643 0.11 0 50.869
## depthError magError type
## 1 99.90 0.12 non-earthquake
## 2 99.90 0.03 non-earthquake
## 3 99.90 0.04 non-earthquake
## 4 99.90 0.02 non-earthquake
## 5 91.03 0.03 non-earthquake
## 6 99.90 0.09 non-earthquake
## 7 99.90 0.06 non-earthquake
## 8 99.90 0.26 non-earthquake
Noto che abbiamo più di 300 osservazioni con valore pari a 99 e sono in gran parte “non-earthquake”. Decidiamo di non rimuoverli per non eliminare dati rilevanti per la classificazione.
## latitude longitude depth mag nst gap dmin rms magNst horizontalError
## 1 19.08667 -155.4027 37.08 2.75 37 195 0.08852 0.18 8 0.89
## depthError magError type
## 1 1.34 4.26 non-earthquake
Decidiamo di rimuovere l’osservazine in questione, dato che la misurazione della magnitudo è molto incerta, migliorando di conseguenza anche il processo di normalizzazione della variabile.
#Outliers più rilevanti rms
#Rms alto corrisponde ad un errore di localizzazione maggiore
earth_noNA |> filter(rms>10)## 95%
## 14.19
## 95%
## 31.61
Notiamo che gli outlier di rms straordinariamente alti corrispondono a valori elevati in horizontalerror e deptherror. Ciò porta quindi ad una bassa precisione dei dati raccolti, quindi optiamo per rimuovere tali osservazioni.
Non è rilevante ai fini dell’analisi, dunque rimuoviamo l’osservazione per permettere una migliore normalizzazione dei dati.
Rimuoviamo i valori dato che una grande dmin potrebbe significare una grande incertezza nelle misurazione, data la difficoltà nella recezione delle onde sismiche.
Nel Dataset abbiamo l’8,4% di Missing Value, possiamo notare inoltre come sembra esserci un pattern ripetuto nei missing value:
Abbiamo molto spesso missing value nelle variabili “gap”, “nst”, “dmin”, “horizonatalError”, “magError” e “magNst”. Queste variabili sono metadati che derivano dalle stime sulla qualità dell’evento sismico, ovvero dal modo in cui sono state calcolate la posizione e la magnitudo dell’evento.
Il fatto che questi campi presentino valori NA (mancanti) quasi sempre insieme può succedere perché il calcolo di questi parametri richiede una copertura minima di stazioni sismiche. Se per un determinato evento non sono disponibili dati sufficienti generalmente anche nelle altre variabili saranno presenti valori mancanti (ad esempio, poche stazioni hanno registrato l’evento o la loro distribuzione non è omogenea).
Visualizziamo come si distribuiscono gli NA
Possiamo notare che, come detto in precedenza, abbiamo molti valori mancanti nella varibili gap(13%), nst(13%), dmin(20%), horizonatalError(18%), magError(15%) e magNst(15%)
Andiamo a visualizzare i pattern degli NA, e come si distribuiscono nelle diversi variabili.
## latitude longitude mag type depth rms depthError nst gap magNst
## 55345 1 1 1 1 1 1 1 1 1 1
## 2316 1 1 1 1 1 1 1 1 1 1
## 491 1 1 1 1 1 1 1 1 1 1
## 919 1 1 1 1 1 1 1 1 1 1
## 349 1 1 1 1 1 1 1 1 1 1
## 36 1 1 1 1 1 1 1 1 1 1
## 37 1 1 1 1 1 1 1 1 1 1
## 203 1 1 1 1 1 1 1 1 1 0
## 12 1 1 1 1 1 1 1 1 1 0
## 2 1 1 1 1 1 1 1 1 1 0
## 1 1 1 1 1 1 1 1 1 0 1
## 1 1 1 1 1 1 1 1 1 0 0
## 1789 1 1 1 1 1 1 1 0 1 1
## 134 1 1 1 1 1 1 1 0 1 1
## 68 1 1 1 1 1 1 1 0 1 0
## 9 1 1 1 1 1 1 1 0 1 0
## 7 1 1 1 1 1 1 1 0 0 1
## 1 1 1 1 1 1 1 1 0 0 1
## 1 1 1 1 1 1 1 1 0 0 1
## 3 1 1 1 1 1 1 1 0 0 1
## 1 1 1 1 1 1 1 1 0 0 0
## 3841 1 1 1 1 1 1 1 0 0 0
## 43 1 1 1 1 1 1 0 1 1 1
## 65 1 1 1 1 1 1 0 1 1 1
## 299 1 1 1 1 1 1 0 1 1 1
## 900 1 1 1 1 1 1 0 1 1 1
## 1710 1 1 1 1 1 1 0 1 1 0
## 2 1 1 1 1 1 1 0 1 0 1
## 2150 1 1 1 1 1 1 0 1 0 0
## 235 1 1 1 1 1 1 0 0 0 1
## 89 1 1 1 1 1 1 0 0 0 0
## 1 1 1 1 1 1 0 1 1 0 1
## 4 1 1 1 1 1 0 1 0 0 1
## 2 1 1 1 1 1 0 1 0 0 1
## 2747 1 1 1 1 1 0 1 0 0 0
## 2 1 1 1 1 1 0 1 0 0 0
## 1 1 1 1 1 1 0 0 1 1 1
## 2 1 1 1 1 1 0 0 1 1 0
## 2 1 1 1 1 1 0 0 1 0 1
## 1 1 1 1 1 1 0 0 1 0 1
## 7 1 1 1 1 1 0 0 0 0 1
## 188 1 1 1 1 1 0 0 0 0 1
## 399 1 1 1 1 1 0 0 0 0 0
## 173 1 1 1 1 0 0 0 0 0 0
## 0 0 0 0 173 3529 6266 9700 9858 11409
## magError horizontalError dmin
## 55345 1 1 1 0
## 2316 1 1 0 1
## 491 1 0 1 1
## 919 0 1 1 1
## 349 0 1 0 2
## 36 0 0 1 2
## 37 0 0 0 3
## 203 0 1 1 2
## 12 0 0 1 3
## 2 0 0 0 4
## 1 1 1 0 2
## 1 0 0 0 5
## 1789 1 1 1 1
## 134 1 0 1 2
## 68 1 0 1 3
## 9 0 1 1 3
## 7 1 1 0 3
## 1 1 0 0 4
## 1 0 1 0 4
## 3 0 0 0 5
## 1 0 1 0 5
## 3841 0 0 0 6
## 43 1 1 1 1
## 65 1 0 1 2
## 299 0 0 1 3
## 900 0 0 0 4
## 1710 0 0 0 5
## 2 0 0 0 5
## 2150 0 0 0 6
## 235 0 0 0 6
## 89 0 0 0 7
## 1 1 1 0 3
## 4 1 1 0 4
## 2 0 0 0 6
## 2747 1 0 0 6
## 2 0 0 0 7
## 1 1 0 1 3
## 2 0 0 0 6
## 2 0 0 1 5
## 1 0 0 0 6
## 7 1 0 0 6
## 188 0 0 0 7
## 399 0 0 0 8
## 173 0 0 0 9
## 11568 13600 15172 81275
Notiamo come diverse volte sono presenti dei valori mancanti in più variabili. Per Facilitare l’imputazione multipla degli NA eliminiamo i casi in cui si hanno missing value in più di 6 variabili: la percentuale di NA non varia molto.
ismiss <- is.na(earth)
ind <- rowSums(ismiss) > 6
indx <- which(ind == T)
earth <- earth[-indx,]
vis_miss(earth, warn_large_data = FALSE)Dividiamo il nostro dataset in train (80%) e test (20%)
#divisione train e test
set.seed(1)
test <- sample(1:nrow(earth), round(0.2*nrow(earth)))
#Train
earth_train <- earth[-test,]
#Classi del Train
class_train <- earth$type[-test]
#Test
earth_test <- earth[test,]
#Classi del test
class_test_NA <- as.factor(earth_test$type)
earth_test <- earth_test |> dplyr::select(-type)Creiamo una funzione per normalizzare le variabili in modo da avere stime più accurate. Calcoliamo il minimo e il massimo per ogni variabile nel training set e trasformiamo il test set utilizzando gli stessi valori.
#Funzione per Normalizzare
normalize_data <- function(train, test) {
numeric_columns <- sapply(train, is.numeric)
train_normalized <- train
test_normalized <- test
train_numeric <- train[numeric_columns]
test_numeric <- test[numeric_columns]
min_vals <- apply(train_numeric, 2, min, na.rm = TRUE)
max_vals <- apply(train_numeric, 2, max, na.rm = TRUE)
# Normalizza i dati di training
train_normalized[numeric_columns] <- sweep(train_numeric, 2, min_vals, FUN = "-")
train_normalized[numeric_columns] <- sweep(train_normalized[numeric_columns], 2, max_vals - min_vals, FUN = "/")
# Normalizza i dati di test usando gli stessi parametri
test_normalized[numeric_columns] <- sweep(test_numeric, 2, min_vals, FUN = "-")
test_normalized[numeric_columns] <- sweep(test_normalized[numeric_columns], 2, max_vals - min_vals, FUN = "/")
return(list(train = train_normalized, test = test_normalized))
}Poichè disponiamo di sufficenti dati nel test set per valutare i nostri modelli procediamo con una strategia passiva per trattare gli NA, ovvero con una casewise delection. In altro modo eseguiamo una strategia attiva per il trattamento dei missing values nel trainig set. Visualizziamo la percentuale di valori mancanti nel training set, che rimane praticamente invariata da quella calcolata in precedenza.
earth_test <- earth[test,]
earth_test <- na.omit(earth_test)
class_test <- as.factor(earth_test$type)
earth_test <- earth_test |> dplyr::select(-type)
vis_miss(earth_train)L’imputazione con PMM è una tecnica semi-parametrica usata per imputare valori mancanti mantenendo la distribuzione originale dei dati. Si inizia stimando un modello di regressione sui dati completi. Per ogni osservazione con un valore mancante, si predice il valore mancante usando il modello. Si trovano poi le osservazioni reali (donatrici) con valori predetti simili e si sceglie casualmente uno dei valori osservati corrispondenti.
#Imputiamo gli NA tramite pmm
my_pred_matrix1 <- 1 - diag(nrow = dim(earth_train)[2], ncol = dim(earth_train)[2])
my_pred_matrix1[,dim(earth_train)[2]] <- 0
earth_train_mice <- mice(earth_train, method = "pmm", seed = 1, printFlag = FALSE, predictorMatrix = my_pred_matrix1)
#Controlliamo come cambiano le distribuzioni con i valori imputati
densityplot(earth_train_mice)Osserviamo come variano le imputazioni multiple delle diverse variabili nei 5 dataset creati. In generale le imputazioni sembrano tutte preservare la distribuzione originale dei dati osservati e inoltre le imputazioni nei 5 dataset sembrano molto simili tra di loro.
Regressione Logistica (pooling) : Utilizziamo glm.mids() che è una versione di glm() che lavora direttamente su oggetti di classe mids, cioè oggetti creati da mice() dopo aver imputato i dati mancanti. Il comando restituisce una lista di modelli glm applicati su ciascun dataset imputato. E poi andremo a valutare tutti i dataset con il comando pool.
#Modello logistico senza imputare gli NA
model <- glm.mids(type ~ depth + mag + nst + gap + dmin + rms + magNst + horizontalError + depthError + magError , data = earth_train_mice, family = "binomial")
pooled <- pool(model)
summary(pooled)Notiamo che praticamente tutti i coefficenti sono significativi, a parte la variabile depthError.
Procediamo con la previsione e la valutazione del modello.
#Previsione con reg. logistica
pred <- lapply(getfit(model), predict, se.fit = T, newdata = earth_test, type = "response")
single_prediction <- sapply(pred, `[[`,"fit")
final_pred <- apply(single_prediction, 1, mean)
final_pred_class <- as.factor(ifelse(final_pred > 0.5, 1, 0))
levels(final_pred_class) <- levels(class_test)
confusionMatrix(final_pred_class, class_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction earthquake non-earthquake
## earthquake 10139 95
## non-earthquake 197 617
##
## Accuracy : 0.9736
## 95% CI : (0.9704, 0.9765)
## No Information Rate : 0.9356
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.7945
##
## Mcnemar's Test P-Value : 3.409e-09
##
## Sensitivity : 0.9809
## Specificity : 0.8666
## Pos Pred Value : 0.9907
## Neg Pred Value : 0.7580
## Prevalence : 0.9356
## Detection Rate : 0.9177
## Detection Prevalence : 0.9263
## Balanced Accuracy : 0.9238
##
## 'Positive' Class : earthquake
##
La previsione nel test porta a buoni risultati come si può vedere dall’Accuracy di 0.9736. Andando nel dettaglio della confusion matrix otteniamo una sensitivity di 0.9809, quindi il modello riesce a predirre molto bene la classe di default (earthquake), data dalla formula della sensitivity.
\[ \text{Sensitivity} = \frac{TP}{TP + FN} \]
D’altro canto la specificity porta un valore più basso, pari a 0.8666. Questo è dovuto probabilmente al forte sbilanciamento delle due classi.
\[ \text{Specificity} = \frac{TN}{TN + FP} \]
Osserviamo che il Pos predicted Value è pari a 0.99057, in accordo con quanto ottenuto precedentemente con la sensitivity.
\[ \text{Pos predicted Value} = \frac{TP}{TP + FP} \]
Osserviamo che il Neg predicted Value è pari a 0.7580, in accordo con quanto ottenuto precedentemente con la specificity.
\[ \text{Neg predicted Value} = \frac{TN}{TN + FN} \]
Il modello predice molto bene la classe earthquake ma non in maniera ottimale la classe più sbilanciata (non-earthquake)
Poichè le classi sono molto sbilanciate utilizziamo la metrica dell’F1-score:
\[ \text{F1-score} = \frac{2}{\frac{1}{P} + \frac{1}{R}} \]
# Funzione per calcolare l'F1-score
f1_score <- function(conf_matrix, positive_class = NULL) {
cm_table <- conf_matrix$table
if (is.null(positive_class)) {
positive_class <- rownames(cm_table)[1] #
}
TP <- cm_table[positive_class, positive_class]
FP <- sum(cm_table[, positive_class]) - TP
FN <- sum(cm_table[positive_class, ]) - TP
precision <- if ((TP + FP) > 0) TP / (TP + FP) else 0
recall <- if ((TP + FN) > 0) TP / (TP + FN) else 0
f1 <- if ((precision + recall) > 0) {
2 * precision * recall / (precision + recall)
} else {
0
}
return(f1)
}
cm <- confusionMatrix(final_pred_class, class_test)
f1_score(cm, positive_class = "non-earthquake")## [1] 0.8086501
Il valore dell’F1-score è di 0.809, questo risultato è collegato a quanto detto in precedenza: questo modello non riesce molto a predirre la classe “non-earthquake” (la classe sbilanciata).
Scegliamo un dataset per andare alla ricerca del K ottimale.
norm <- normalize_data(complete(earth_train_mice,1),earth_test)
train_std <- norm$train
test_std <- norm$testDividiamo il training set in validation e training per valutare la migliore accuracy per i diversi valori di k (con k da 2 a 15)
train_std <- train_std |> dplyr::select(-c("type","latitude","longitude"))
set.seed(1)
validation <- sample(1:nrow(train_std),0.2*nrow(train_std))
validation_set <- train_std[validation,]
class_validation <- class_train[validation]
earth_train_complete2 <- train_std[-validation,]
class_train2 <- class_train[-validation]
calc_class_err = function(actual, predicted) { mean(actual != predicted) }
k <- c(2:15)
accuracy <- vector()
err_k = rep(x = 0, times = length(k))
result <- matrix(NA, ncol = length(k), nrow = 2)
for(i in 1:length(k)){
pred <- as.factor(knn(earth_train_complete2,validation_set,class_train2, k = k[i]))
accuracy[i] <- mean(pred == class_validation)
result[1,i] <- k[i]
result[2,i] <- accuracy[i]
err_k[i] = calc_class_err(class_validation, pred)
}
k[which.max(result[2,])] ## [1] 5
Otteniamo un k ottimale uguale a 5.
plot(k,result[2,] , type = "b", col = "dodgerblue", cex = 1, pch = 20,
xlab = "k, number of neighbors", ylab = "Accuracy",
main = "Accuracy vs Neighbors")Utilizziamo la seguente strategia di pooling: stimiamo un modello KNN per ogni dataset creato dell’imputazione e teniamo come predizione finale la classe più predetta per ogni osservazione nei 5 modelli.
predizioni <- matrix(NA, nrow = nrow(earth_test), ncol = 5)
for(i in 1:5){
std <- normalize_data(complete(earth_train_mice,i) |> dplyr::select(-c("type","latitude","longitude")), earth_test|> dplyr::select(-c("latitude","longitude")))
train <- std$train
test <- std$test
pred <- as.factor(knn(train, test = test, class_train, k = 3))
predizioni[,i] <- pred
}
elemento_piu_frequente <- function(riga){
tabella <- table(riga)
moda <- names(tabella)[tabella == max(tabella)]
return(moda[1]) # in caso di parità, prende il primo
}
pool_prediction <- apply(predizioni, 1, elemento_piu_frequente)
pool_prediction <- as.factor(pool_prediction)
levels(pool_prediction) <- levels(class_test)
confusionMatrix(pool_prediction,class_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction earthquake non-earthquake
## earthquake 10198 112
## non-earthquake 138 600
##
## Accuracy : 0.9774
## 95% CI : (0.9744, 0.9801)
## No Information Rate : 0.9356
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8155
##
## Mcnemar's Test P-Value : 0.1138
##
## Sensitivity : 0.9866
## Specificity : 0.8427
## Pos Pred Value : 0.9891
## Neg Pred Value : 0.8130
## Prevalence : 0.9356
## Detection Rate : 0.9231
## Detection Prevalence : 0.9332
## Balanced Accuracy : 0.9147
##
## 'Positive' Class : earthquake
##
La previsione nel test porta a buoni risultati come si può vedere dall’Accuracy di 0.9774. Il modello sembra migliorare leggermente rispetto alla regressione logistica. Andando nel dettaglio della confusion matrix otteniamo una sensitivity di 0.9863, quindi il modello riesce a predirre molto bene la classe di default (earthquake), data dalla formula della sensitivity.
\[ \text{Sensitivity} = \frac{TP}{TP + FN} \]
D’altro canto la specificity porta un valore più basso, pari a 0.8427. Il modello peggiora leggermente da questo punto di vista rispetto alla regressione logistica.
\[ \text{Specificity} = \frac{TN}{TN + FP} \]
Osserviamo che il Pos predicted Value è pari a 0.9866, in accordo con quanto ottenuto precedentemente con la sensitivity.
\[ \text{Pos predicted Value} = \frac{TP}{TP + FP} \]
Osserviamo che il Neg predicted Value è pari a 0.8130, in accordo con quanto ottenuto precedentemente con la specificity. Il knn quindi migliora molto il Negative Pred Value rispetto alla regressione logistica.
\[ \text{Neg predicted Value} = \frac{TN}{TN + FN} \]
## [1] 0.8275862
In questo caso il valore dell’F1_score migliora leggermente rispetto alla regressione logistica, il modello KNN riesce a predirre meglio la cllasse sblilanciata ma comunque non in modo ottimale.
Nonostante le assunzioni della QDA non sembrano verificate dagli istogrammi visualizzati in precedenza andiamo a stimare lo stesso il modello.
pred_qda1 <- predict(earth_qda1, na.omit(earth_test))
post_qda <- apply(pred_qda1$posterior, 1, max)
confusionMatrix(pred_qda1$class , class_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction earthquake non-earthquake
## earthquake 8917 19
## non-earthquake 1419 693
##
## Accuracy : 0.8698
## 95% CI : (0.8634, 0.8761)
## No Information Rate : 0.9356
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.4365
##
## Mcnemar's Test P-Value : <2e-16
##
## Sensitivity : 0.8627
## Specificity : 0.9733
## Pos Pred Value : 0.9979
## Neg Pred Value : 0.3281
## Prevalence : 0.9356
## Detection Rate : 0.8071
## Detection Prevalence : 0.8088
## Balanced Accuracy : 0.9180
##
## 'Positive' Class : earthquake
##
Si può vedere subito come la QDA peggiora notevolmente rispetto agli altri modelli, questo poteva essere anticipato dalla verifica delle assunzioni della QDA.
###Previsioni con dati senza NA
earth_lda_tst_pred1 <- predict(earth_lda, earth_test)
post_lda <- apply(earth_lda_tst_pred1$posterior, 1, max)
confusionMatrix(earth_lda_tst_pred1$class,class_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction earthquake non-earthquake
## earthquake 10038 526
## non-earthquake 298 186
##
## Accuracy : 0.9254
## 95% CI : (0.9204, 0.9302)
## No Information Rate : 0.9356
## P-Value [Acc > NIR] : 1
##
## Kappa : 0.2731
##
## Mcnemar's Test P-Value : 2.617e-15
##
## Sensitivity : 0.9712
## Specificity : 0.2612
## Pos Pred Value : 0.9502
## Neg Pred Value : 0.3843
## Prevalence : 0.9356
## Detection Rate : 0.9086
## Detection Prevalence : 0.9562
## Balanced Accuracy : 0.6162
##
## 'Positive' Class : earthquake
##
Possiamo notare un accuracy alta, causata dallo sbilanciamento delle classi. infatti se ci focalizziamo su indici più specifici notiamo che la specificity è pari a 0.2612 e il negative predicted value è pari a 0.3843. Analogamente alla QDA questi risultati erano prevedibili dalla verifica delle assunzioni.
tree.type <- tree(type ~ depth + mag + nst + gap + dmin + rms + magNst + horizontalError + depthError + magError, earth_train)
plot(tree.type,type = "uniform")
text(tree.type, ,pretty = 0,cex = 0.6)Andiamo adesso a stimare la miglior grandezza dell’albero tramite la funzione cv.tree().
#CV dell'alfa
set.seed(1)
cv.type <- cv.tree(tree.type, FUN = prune.misclass)
ggplot(data = data.frame(size = cv.type$size, dev = cv.type$dev)) +
geom_line(aes(x = size, y = dev), col = "#31A354") + geom_point(aes(x = size, y = dev),alpha = 0.2, size = 4, col = "#31A354") +
geom_hline(yintercept = min(cv.type$dev), col = "#8C6D31", linetype = "dashed") Notiamo come le dimensioni 7 e 8 dell’albero ottengono lo stesso risultato che minimizzano la devianza.
Andiamo quindi a ristimare l’albero con la grandezza 7 anzichè 8, scegliamo la dimensione minore perchè in un nodo si ottiene la stessa risultante, come si può vedere nel grafico precedentemente riportato.
best_k <- 7
pruned.tree <- prune.tree(tree.type, best = best_k)
plot(pruned.tree,type = "uniform")
text(pruned.tree, ,pretty = 0,cex = 0.6)#Previsione con decision trees
tree.pred <- predict(pruned.tree , earth_test , type = "class")
confusionMatrix(tree.pred, class_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction earthquake non-earthquake
## earthquake 10247 78
## non-earthquake 89 634
##
## Accuracy : 0.9849
## 95% CI : (0.9824, 0.9871)
## No Information Rate : 0.9356
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.8755
##
## Mcnemar's Test P-Value : 0.439
##
## Sensitivity : 0.9914
## Specificity : 0.8904
## Pos Pred Value : 0.9924
## Neg Pred Value : 0.8769
## Prevalence : 0.9356
## Detection Rate : 0.9275
## Detection Prevalence : 0.9346
## Balanced Accuracy : 0.9409
##
## 'Positive' Class : earthquake
##
La previsione nel test porta a buoni risultati come si può vedere dall’Accuracy di 0.9849. Andando nel dettaglio della confusion matrix otteniamo una sensitivity di 0.9914, quindi il modello riesce a predirre molto bene la classe di default (earthquake), migliore degli altri modelli
D’altro canto la specificity porta un valore più basso, pari a 0.8904, acnhe in questo caso i risultati ottenuti sono migliori degli altri modelli.
Osserviamo che il Pos predicted Value è pari a 0.9924, in accordo con quanto ottenuto precedentemente con la sensitivity.
Osserviamo che il Neg predicted Value è pari a 0.8769, anche questo valore è migliorato di parecchio in confronto agli altri modelli stimati precedentemente.
Valutiamo congiuntamente tutti i risultati dei modelli applicati attraverso la Curva di ROC. La curva ROC è un grafico che mette in relazione la sensibilità (True Positive Rate) con la specificità (1 − False Positive Rate) per diverse soglie di classificazione. Serve per valutare quanto bene il modello distingue tra le due classi. Più la curva è vicina all’angolo in alto a sinistra, migliore è il modello e maggiore sara l’area sotto la curva (AUC), la quale quantifica la qualità della classificazione.
#Salvo le posterior del pooling del KNN
posterior <- matrix(NA, nrow = nrow(earth_test), ncol = 5)
for(i in 1:5){
std <- normalize_data(complete(earth_train_mice,i) |> dplyr::select(-type), earth_test)
train <- std$train
test <- std$test
pred <- knn(train, test = test, class_train, k = 3, prob = T)
prob_attr <- attr(pred, "prob")
posterior[,i] <- prob_attr
}
posterior_pool <- rowMeans(posterior)
prob_Earthquake <- ifelse(pool_prediction == "earthquake", posterior_pool, 1 - posterior_pool)Calcoliamo la curva di ROC per i vari modelli e la rispettiva area sottostante la curva (AUC).
#per il knn
roc_obj <- roc(response = class_test, predictor = prob_Earthquake, levels = c("earthquake","non-earthquake"))## Setting direction: controls > cases
#per la logistica
roc_log <- roc(response = class_test, predictor = final_pred, levels = c("earthquake","non-earthquake"))## Setting direction: controls < cases
roc_albero <- roc(response = class_test, predictor = post_tree, levels = c("earthquake","non-earthquake"))## Setting direction: controls > cases
roc_qda <- roc(response = class_test, predictor = post_qda, levels = c("earthquake","non-earthquake"))## Setting direction: controls > cases
roc_lda <- roc(response = class_test, predictor = post_lda, levels = c("earthquake","non-earthquake"))## Setting direction: controls > cases
# Estrazione dei punti per il grafico
roc_df <- data.frame(
fpr = rev(1 - roc_obj$specificities), # false positive rate = 1 - specificity
tpr = rev(roc_obj$sensitivities) # true positive rate = sensitivity
)
df_log <- data.frame(
fpr = rev(1 - roc_log$specificities),
tpr = rev(roc_log$sensitivities)
)
df_albero <- data.frame(
fpr = rev(1 - roc_albero$specificities),
tpr = rev(roc_albero$sensitivities)
)
df_qda <- data.frame(
fpr = rev(1 - roc_qda$specificities),
tpr = rev(roc_qda$sensitivities)
)
df_lda <- data.frame(
fpr = rev(1 - roc_lda$specificities),
tpr = rev(roc_lda$sensitivities)
)
auc_knn <- auc(roc_obj)
auc_log <- auc(roc_log)
auc_albero <- auc(roc_albero)
auc_qda <- auc(roc_qda)
auc_lda <- auc(roc_lda)
# curva interattiva con plotly
plot_ly() %>%
add_trace(
data = roc_df,
x = ~fpr,
y = ~tpr,
type = 'scatter',
mode = 'lines',
name = paste0("kNN (AUC = ", round(auc_knn, 3), ")"),
line = list(color = 'blue'),
hoverinfo = "text+name"
) %>%
add_trace(
data = df_log,
x = ~fpr,
y = ~tpr,
type = 'scatter',
mode = 'lines',
name = paste0("Logistica (AUC = ", round(auc_log, 3), ")"),
line = list(color = 'red'),
hoverinfo = "text+name"
) %>%
add_trace(
data = df_albero,
x = ~fpr,
y = ~tpr,
type = 'scatter',
mode = 'lines',
name = paste0("Albero Decisionale (AUC = " , round(auc_albero, 3), ")"),
line = list(color = 'green'),
hoverinfo = "text+name"
) %>%
add_trace(
data = df_qda,
x = ~fpr,
y = ~tpr,
type = 'scatter',
mode = 'lines',
name = paste0("QDA (AUC = ", round(auc_qda, 3), ")"),
line = list(color = 'purple'),
hoverinfo = "text+name"
) %>%
add_trace(
data = df_lda,
x = ~fpr,
y = ~tpr,
type = 'scatter',
mode = 'lines',
name = paste0("LDA (AUC = ", round(auc_lda, 3), ")"),
line = list(color = 'yellow'),
hoverinfo = "text+name"
) %>%
layout(
title = "Confronto Curve di ROC",
xaxis = list(title = "False Positive Rate"),
yaxis = list(title = "True Positive Rate"),
legend = list(x = 0.8, y = 0.2)
)Osservando il grafico e il valore dell’AUC, i modelli migliori risultano essere il KNN con AUC = 0.948 e la Regressione logistica con un AUC = 0.981. L’albero decisionale presenta comunqe un AUC pari a 0.893.
In conclusione dal punto di vista della curva di ROC il modello migliore sembra essere la Regressione Logistica, ma osservando la confusion matrix notiamo la sua scarsa capacità di identificare la classe meno rappresentata (non-earthquake). Il discorso è analogo per il KNN, mentre l’albero decisionale nonostante abbia un valore inferiore nella curva di ROC ottiene risultati migliori nella confusion matrix (specificity, neg pred value).